Rotating Spiral Waves with Phase-Randomized Core in Non-locally Coupled 

Oscillators 



m 
O 
O 

(N 
Or 

OS 



c/3 

PL, 

d 



> 

(N 

in 
o 
o> 
o 
m 
o 



X 



Shin-ichiro Shim£0 and Yoshiki Kuramoto 
Department of Physics, Graduate School of Sciences, Kyoto University, Kyoto 606-8502, JAPAN 

(Dated: February 8, 2008) 

Rotating spiral waves with a central core composed of phase-randomized oscillators can arise 
in reaction-diffusion systems if some of the chemical components involved are diffusion- free. This 
peculiar phenomenon is demonstrated for a paradigmatic three-component reaction-diffusion model. 
The origin of this anomalous spiral dynamics is the effective non-locality in coupling, whose effect 
is stronger for weaker coupling. There exists a critical coupling strength which is estimated from a 
simple argument. Detailed mathematical and numerical analyses are carried out in the extreme case 
of weak coupling for which the phase reduction method is applicable. Under the assumption that the 
mean field pattern keeps to rotate steadily as a result of a statistical cancellation of the incoherence, 
we derive a functional self-consistency equation to be satisfied by this space-time dependent quantity. 
Its solution and the resulting effective frequencies of the individual oscillators are found to agree 
excellently with the numerical simulation. 

PACS numbers: 82.40.Ck, 05.45.Xt 



I. INTRODUCTION 



Rotating spiral waves represent a most universal form 
of patterns appearing in reaction-diffusion systems and 
other dissipative media of oscillatory and excitable na- 
ture. Most of the recent experimental and theoretical 
studies on rotating spiral waves in reaction-diffusion sys- 
tems have focused on their complex behavior such as 
core meandering in 2D0, 0], manipulation of the pat- 
tern using photo-sensitive BZ reaction^ 3, and the 
topology and dynamics of the singular filaments of 3D 
scroll waves0, El 0- Slightly deviated from this main- 
stream, the possibility of a new type of spiral dynam- 
ics caused by a universal mechanism was proposed re- 
cently by the present authors 8J. This is characterized 
by the appearance of a local group of oscillators near 
the center of spiral rotation where the oscillators be- 
have individually rather than collectively. A simple class 
of three-component reaction-diffusion systems involving 
two diffusion-free components and an extra diffusive com- 
ponent proved to exhibit this type of anomaly. Here the 
last component plays the role of a coupling agent allow- 
ing the otherwise independent local oscillators to commu- 
nicate with each other, where the communication takes 
place non-locally. The crucial parameter to this peculiar 
spiral dynamics is the strength of the non-local coupling. 
If it is sufficiently large, the characteristic wavelength of 
the pattern, especially the radius of the spiral core, be- 
comes longer than the coupling radius. Consequently, 
the coupling becomes effectively local, i.e., diffusive, and 
there is nothing peculiar about the resulting spiral pat- 
tern. As the coupling becomes weaker, in contrast, its 
non-local nature becomes stronger, and finally a small 
group of phase-randomized oscillators starts to be created 
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near the center of rotation. We find in the present pa- 
per that under certain conditions the phase-randomized 
core is stationary in a statistical sense. This allows us to 
formulate a statistical theory with which the entire sys- 
tem dynamics, collective and individual, can completely 
be specified. Actually, the exact self-consistent theory 
developed here provides a rare example of statistical the- 
ories associated with large systems of limit-cycle oscilla- 
tors when spatial degrees of freedom are involved. 

The organization of the present paper is the follow- 
ing. In Sec. II, we start with a brief review of some gen- 
eral features of the three-component reaction-diffusion 
model introduced earlier, and show how it is reduced to a 
two-component system of non-locally coupled oscillators. 
Then, adopting a specific model for the local oscillators, 
we present some results of our numerical simulation re- 
vealing the fact that the spiral core can be coherent or 
incoherent depending on the coupling parameter. We 
shall also see that the critical coupling strength associ- 
ated with the onset of incoherence can be estimated from 
a simple argument. Section III and IV, each devoted 
to numerical and mathematical analyses, are concerned 
with the special situation where the coupling is suffi- 
ciently weak. Then the so-called phase reduction method 
is applicable, by which a phase oscillator model with non- 
local coupling is derived. What is remarkable is the fact 
that, as opposed to the conventional view, description 
of the spiral dynamics in terms of the phase oscillator 
model does not lead to a topological contradiction but 
can even provide its precise description. Using this non- 
local phase model, we develop a mean field theory simi- 
lar to Kuramoto's f 975 theory on the onset of collective 
synchronization in globally coupled oscillators 9J . The 
present mean field theory is based on the assumption of 
steady rotation of the mean field pattern. Owing to this 
assumption, we can derive a functional self-consistency 
equation to be satisfied by the mean field. Numerical 
solution of this functional equation is confirmed to agree 
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exceedingly well with the simulation results. 

II. REACTION-DIFFUSION MODEL AND ITS 
REDUCED FORM 

A. Effective Non-locality in Reaction-Diffusion 

Systems 

Consider a three-component reaction-diffusion system 
of the following form. 

d t X(r,t) = f(X,Y) + K(B-X), (1) 
d t Y(r,t)=g(X,Y), (2) 
Td t B{r,t) = -B + DV 2 B + X. (3) 

The system is supposed to extend sufficiently in two di- 
mensions. The above model has recently been used as 
a paradigmatic model for the study of various aspects 
of self-oscillatory fields where the effective non-locality 
in coupling plays a crucial role [3, El El- The first two 
equations with K = represent a local limit-cycle os- 
cillator. Our system may therefore be interpreted as a 
continuum limit of a large assembly of oscillators with- 
out direct mutual coupling which are suspended in a dif- 
fusive chemical with concentration B. The last quantity 
plays the role of a coupling agent only by which the local 
oscillators can mutually communicate. For simplicity, a 
cross coupling between the local oscillators and the dif- 
fusive field has been introduced in a linear form and only 
between X and B. The coupling term K(B — X) may 
equivalently be replaced with a more natural form KB if 
f(X, Y) is suitably redefined, but we will work with the 
first form for its mathematical convenience to be seen 
later. 

Our system, possibly with various modifications and 
generalizations, bears some resemblance to biological 
populations of oscillatory and excitable cells such as 
suspensions of yeast cells under glycolysis and slime 
mold amoebae in a certain phase of their life cycle 
IT3I Il4| . One may also note some similarity of the 
above model to the recently developed version of the 
Belousov-Zhabotinsky reaction using water-in-oil AOT 
micro-emulsion |l5l lift Il7j. Some of the interesting the- 
oretical aspects of our reaction-diffusion model have al- 
ready been reported H, El E] • 

When the characteristic time scale of B, denoted by 
r, is sufficiently small, this component can be eliminated 
adiabatically by solving the equation 

= -B + DV 2 B + X. (4) 

The solution of the above equation is expressed in terms 
of the Green function G(r) in the form 

B(r,t) = J G(r-r')X(r',t)dr'. (5) 

If our system is infinitely extended, G(r) is radially sym- 
metric, and for spatial dimension two it is given by a 



modified Bessel function of the second kind with the char- 
acteristic length scale ro = V"D, i.e., 



Note that the above G(r) satisfies the normalization con- 
dition J G(r)dr = 1, and behaves asymptotically as 
G(r) r~j exp(— r/ro)/y / r/ro for r ro. We may call 
B(r,t) the space-time dependent mean field because this 
quantity roughly represents a mean value of X(r' , t) over 
a circular domain with the radius of 0(VT)) centered at 
r. Our original reaction-diffusion system has now been 
reduced to the system of Eqs. QJ, J2J and (JSJ), which rep- 
resents a two-component oscillatory field with non-local 
coupling. 

Suppose that we change parameter K, which measures, 
in terms of the reduced system, the strength of the non- 
local coupling. If K is sufficiently large, the characteris- 
tic wavelength of the pattern, denoted by l p , will be far 
longer than the coupling radius (as is justified below). 
Then the long-wavelength approximation can applied to 
Eq. ©, giving K(B - X) ~ DV 2 X, where D = KD. 
Thus, the non-local coupling practically reduces to a dif- 
fusive coupling in this strong-coupling case. One may 
check the consistency of the above argument by noticing 
the fact that the result of this diffusion-coupling approx- 
imation itself tells that l p estimated from Eq. (JJJ scales 

like l p ~ \JD = V KD. Thus, sufficiently large K implies 
lp larger than the coupling radius \fD, so that our long- 
wavelength assumption proves to be consistent. In any 
case, our system for strong coupling reduces to a two- 
component reaction-diffusion system, which, however, is 
of our little concern in the present paper. 

The situation of our interest is the opposite case in 
which K is so small that l p becomes comparable with or 
even smaller than the coupling radius of 0(>/~D). Then 
the diffusion-coupling approximation breaks down, and 
the system comes to behave in an unusual manner. It 
should be noted here that the evolution equations them- 
selves are free of characteristic length scale below the 
coupling radius. Therefore, once l p comes to fall within 
the coupling radius, or equivalently, once spatial varia- 
tions with wavelengths smaller than the coupling radius 
are generated spontaneously, then there is no reason why 
spatial variations of even smaller wavelengths should not 
occur. We suspect therefore that the kind of anomaly of 
our concern might be characterized by a fragmentation 
of the pattern down to infinitesimal spatial scales. This 
is actually the case, which we show below by presenting 
some numerical results on the non-locally coupled system 
given by Eqs. QJ, © and JSJ with specific forms for / 
and g. 
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B. Case of the FitzHugh-Nagumo Oscillators 

As a simple model for the local oscillators, let us con- 
sider the FitzHugh-Nagumo model given by 

f = a- 1 {(X-X 3 )-Y}, g = aX + b. (7) 

We fix the parameter values as a = 1.0, b = 0.2 and 
a = 0.1, so that the system is well in the self-oscillatory 
regime. 

We carried out a numerical analysis on a non-locally 
coupled field of oscillators described by Eqs. CJ, © ; © 
and l (7| l. The system is defined over the square domain 
x, y G [0, L] where B satisfies the free boundary con- 
ditions, namely, the vertical component of V B to the 
boundaries vanishes. Thus, the Green function G differs 
in this case from the form given by Eq. © especially 
near the boundaries. In practical numerical simulation, 
our continuous space was replaced with a square lattice 
of oscillators of N x N lattice points, a typical value of 
N being 2048. At each time step, B was calculated from 
Eq. ©, or equivalently Eq. (@J, by means of a spatial 
Fourier transform. The 4th order Runge-Kutta scheme 
was adopted for the time-integration of Eqs. © and J2J. 

Some numerical results for two representative values of 
K are illustrated in FIG.Q] Figure^(a) to (c) correspond 
to the case of large K. They respectively show an overall 
spiral pattern, its blowup near the center of rotation, and 
the phase portrait of the pattern in the X-Y plane. The 
last quantity is given by a set of N 2 points in the X-Y 
plane each representing the state of a local oscillator at 
a given time. In usual reaction-diffusion systems such 
as those modeled with two-component reaction-diffusion 
equations, the phase portrait associated with a spiral pat- 
tern is considered to form a simply connected object in- 
volving a phase singularity. This is a natural consequence 
of the homeomorphism which is supposed to characterize 
the mapping between the physical space and the state 
space. The same property seems to hold in the present 
case of large K, and this is consistent with the fact al- 
ready noted that for sufficiently strong coupling our sys- 
tem reduces to a two-component reaction-diffusion sys- 
tem. 

Figure ^(d) to (f) correspond to the case of small K. 
The overall spiral pattern does not seem qualitatively 
different from that for large K. As is clear from Fig. 
(e), however, closer observation of the core structure re- 
veals a completely new feature of the pattern. This is the 
appearance of a group of oscillators near the center of ro- 
tation where the oscillators seem to behave individually 
rather than collectively. The corresponding phase por- 
trait, which is presented in FIG. H(f), no longer seems 
to tend to a simply connected object in the continuum 
limit. The hole created in the phase portrait gives a clear 
indication of the breakdown of the homeomorphism men- 
tioned above. It may alternatively be said that a pair of 
local oscillators situated infinitely close to each other are 
not always so close in the state space, which says noth- 
ing but a loss of spatial continuity of the pattern. At the 



same time, the phase singularity, which is generally con- 
sidered as a central characteristic shared by spiral pat- 
terns, seems to be lost, i.e., the pattern no longer seems 
to involve a special local oscillator for which the phase 
cannot be defined. 

The origin of the spiral core anomaly of this kind may 
qualitatively be understood in the following way. Our 
primary question is why the core region is the most frag- 
ile part of the pattern with respect to the collapse of 
spatial continuity. In order to see why, it is convenient to 
look upon Eqs. © and J2Jl as describing a single oscillator 
driven by a forcing field B whose spatial variation is ex- 
pected to be relatively smooth from its definition given by 
Eq. (0} . Wherever the oscillation amplitude of B is suffi- 
ciently large, the oscillators will individually synchronize 
with the motion of B, so that a local group of such oscil- 
lators will mutually synchronize also. The corresponding 
local pattern will then look continuous and smooth. This 
is considered to be the case for those oscillators far apart 
from the central core, because the oscillation amplitude 
of B there should be relatively large. In contrast, close 
to the central part of the pattern, where the oscillation 
amplitude of B should be relatively small, synchroniza- 
tion becomes more difficult. Loss of mutual synchrony 
implies the appearance of a group of phase-randomized 
oscillators. 



C. Estimation of K c 

From the numerical data presented in FIG.^ one may 
expect the existence of a critical value of the coupling 
strength, denoted as K c , associated with the onset of in- 
coherence. We now try to estimate K c for our system of 
non-locally coupled FitzHugh-Nagumo oscillators given 
by Eqs. (QJ , J5J , © and 0, where the spatial exten- 
sion is supposed to be infinite. Consider first the situ- 
ation where the coupling is large enough for the system 
to sustain a rigidly rotating spiral wave with sufficient 
spatial smoothness. The corresponding solution is repre- 
sented by A s (r,t) = (X s (r,t),Y s (r,t)). Let the center 
of rotation be at r = 0. By assumption, the oscillator 
right at the center is motionless, i.e., A s (0, t) — (X c , Y c ), 
where X c and Y c are time- independent. Our question is 
at which value of K this fixed point becomes unstable 
and the oscillator there starts to oscillate. To consider 
this problem, it is convenient to work with the afore- 
mentioned mean-field picture by which we look upon the 
local oscillators as being subject to a common space-time 
dependent field B. The mean field pattern should also 
rotate rigidly around r = 0, so that the central oscilla- 
tor is subject to a constant forcing B(0). The system 
of Eqs. Q and @ describing this particular oscillator 
form an autonomous two-dimensional dynamical system, 
so that once B(0) is known the value of the fixed point 
(X c , Y c ) and its stability will easily be found. The value 
of B(0) can actually be estimated from Eq. JSJ by devel- 
oping X s (r,t) into a Taylor series about r = 0, which is 
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FIG. 1: Spiral patterns exhibited by non-locally coupled 
FitzHugh-Nagumo oscillators for two representative cases of 
strong coupling ((a),(b) and (c), K — 10.0) and weak cou- 
pling ((d), (e) and (f), K = 5.0). Other parameters are fixed 
as a = 1.0, b = 0.2, a = 0.1 and D = 1. (a) and (d): Overall 
patterns of the component X displayed in gray scale, (b) and 
(e): Their structures near the core, (c) and (f): Correspond- 
ing phase portraits in the X-Y plane, where the nullclines 
f(X, Y) — and g(X, Y) — are also indicated. 



allowed owing to the assumed smoothness of the pattern. 
It is clear that, as a result of the isotropy of the cou- 
pling function G, there is no contribution to B(0) from 
the first order expansion terms. If the contribution from 
the second order terms is negligible, i.e., if the nonlinear 
variation of X s within the coupling range about r = 
is negligible, then we may simply put B(0) — X c . With 
this approximation, it is clear from Eqs. and J2J that 
the fixed point (X c , Y c ) is identical with the intersection 
of the nullclines / = g = 0, i.e., the unstable fixed point 
of the local oscillators. Its linear stability is also easy to 
analyze. The result is that the critical coupling strength 
is given by K c = (1 — 3X,?) ja below which the fixed point 
(X C ,Y C ) becomes Hopf unstable. Applying the values of 
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FIG. 2: Spiral patterns exhibited by non-locally coupled 
FitzHugh-Nagumo oscillators with K = 2.0, presented in a 
similar manner to FIG.0 For this value of K, the amplitude 
degrees of freedom become almost irrelevant. 

a, b and a used in our numerical simulations, we obtain 
K c = 8.8. This value of K c is consistent with our di- 
rect numerical simulation, although its precise numerical 
determination is yet unavailable. 

III. SPIRAL DYNAMICS IN NON-LOCALLY 
COUPLED PHASE OSCILLATORS 

In order to look into the nature and origin of our 
anomalous spiral dynamics in further detail, we now con- 
sider the situation where the coupling is much weaker 
than K c . Figure [5] shows some results from our numer- 
ical simulation carried out for K = 2.0, presented in a 
similar manner to FIG.Q] While there seems nothing un- 
usual about the overall spiral pattern, the corresponding 
phase portrait forms a ring with a relatively thin periph- 
ery, which is totally unlike a simply connected object. 
We may alternatively say that the oscillator amplitude 
everywhere takes almost a full value. This is exactly the 
situation where the so-called phase description is applica- 
ble. In fact, as we see later in this section, a simple phase 
oscillator model with non-local coupling can develop a 
spiral pattern with phase-randomized core similar to the 
above. 

We now present a brief review of the phase reduc- 
tion method m the form appropriate for the present 
purposes. Each of our local oscillators without cou- 
pling is described by a two-dimensional dynamical sys- 
tem dA/dt = F(A), where A = (X, Y) and F = (f,g). 
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Let its stable time-periodic solution with frequency w 
be given by A (u>t) — (X (ujt) 7 Y (ujt)), which is a 2n- 
periodic function of uit. The corresponding limit-cycle 
orbit is represented by C. Phase tj> associated with this 
oscillator must be defined outside C as well as on C. Most 
conveniently, it is defined in such a way that the free mo- 
tion of the oscillator satisfies d<p/dt = ui regardless of 
initial conditions. This requires that as a scalar field 
4>{A) satisfies the identity grad^-i^A) = ui. The whole 
X-Y plane is then filled with equi-phase lines which are 
called isochrons, one of which is chosen to correspond to 
the zero phase. Corresponding to each Rvalue, a point 
Ao((f>) on C is determined uniquely which says nothing 
but the fact that an isochron and C intersect at a single 
point. 

When the non-local coupling is introduced, the equa- 
tion for each local oscillator is modified as 

d t A(r,t)=F(A)+p(r,t), (8) 

where 




FIG. 3: Phase-coupling function T(<j>) versus <j> for coupled 
FitzHugh-Nagumo oscillators. This quantity can be used for 
the study of non-locally coupled phase oscillators given by 
Eq. JTUJ. 



p(r,t) = (p x (r,t),0), 

px{r,t) = K J G(r-r')[X(r',t)-X(r,t)\ dr'. 

Correspondingly, the equation for the phase is modified 
as 

d t (j) = gr&& A 4> ■ {F{A) + p) = uj + grad A • p. (9) 

If the perturbation p is sufficiently weak, which we as- 
sume now, the oscillator will keep staying on C in good 
approximation. Then grad^c^ in Eq. @ may safely be 
evaluated on C, or 



where 



grad A </>~(Z x (</,),Z y (<«), 



Zx{4>) = [dx<p(A)] A=AoW , 



and Zy{4>) is defined similarly. At the same time, px 
may be approximated with 

p x ~K f G(r-r')[X (cf>{r',t))-X (<j>{r,t))]dr'. 



Thus, the phase equation becomes 
d t <f>(r,t) = u> + KZ x {4>{r,t)) 

x / G{r-r')[X a {<j>{r',t))-X {<j>{r,t))]dr'. 



Since the small effect of the perturbation on dt4> can be 
time-averaged over one cycle of oscillation 18j, the phase 
equation finally takes the form 



d t (t>(r, t) = lu + K J G(r - r')T(<f>(r, t) - <j>(r' , t))dr', 

(10) 



where 



2tt 



2tt 



</>') = — Z x (\+(f>) [X (X + <f>') - X (X + <P)] dX. 



By using the above formula, T(<fi) may be computed nu- 
merically if the forms of / and g are given explicitly. For 
the present case of FitzHugh-Nagumo oscillators, numer- 
ically obtained r(^) is displayed in FIG. [21 

The phase-coupling function r(0), which is a 2ir- 
periodic function of </>, generally involves various harmon- 
ics, and this is also true of the curve given in FIG. 
We still expect that the spiral dynamics of our concern 
does not depend so heavily on the specific form of r(<^>). 
Therefore, in order for further mathematical analysis to 
be practicable, we will work with a simplest in-phase type 
coupling function, i.e., T(cj)) = — sin(0 + a) (\a\ < tt/2). 
Thus, the phase equation takes the form 

d t ^>{r,t)=u)-K / G(r-r')sm((f>(r ) t)-(f>(r' ) t)+a)dr' 



(11) 

for which an in-depth mathematical analysis is possible 
as we see below. 

Before proceeding to the analysis of Eq. ifTTl) . we 
remark that the above phase equation is also a cor- 
rect reduced form of a non-local version of the complex 
Ginzburg-Landau eauation|19j. the latter itself being a 
reduced form of our three-component reaction-diffusion 
model close to the Hopf bifurcation and comparably close 
to the limit of vanishing coupling^]] . This fact gives 
a further support to our view that the application of 
Eq. to our problem is reasonable. 

We are still far from a full understanding of the solution 
to the universal equation and our concern below is 
its spiral wave solution in two dimensions. Although the 
equation involves four parameters lo, K, r$ and a, the 
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only relevant parameter is a. The reason is the following. 
Firstly, r , on which G(r) depends (see Eq. ©), may 
be chosen to be the length unit, so that we may put 
?"o = 1- Similarly, the coupling strength K may be fixed 
to 1 by suitably choosing the time unit. The natural 
frequency ui can be eliminated by working with a suitable 
co-moving frame of reference, i.e., via the transformation 
<f> — > <f> + ujt. In the following analysis, however, the 
irrelevant parameter u> is retained as a nonzero constant, 
while we choose a = 0.3 and r Q = K = I. 

Numerical simulation of Eq. Ijllfl was carried out in a 
two-dimensional system. The numerical scheme adopted 
is the same as that explained in the previous section. As 
expected, we see from FIG. 0] the appearance of rotating 
spiral waves with a disordered group of oscillators near 
the core very similar to what we have seen in the previous 
section. 

For the arguments developed below, it is convenient to 
define a mean field W(r,t) through 

W{r,t) = J G(r-r')exp[i</)(r',t)]dr'. (12) 

The modulus R and the phase of this complex quantity 
are defined by 

W{r, t) = R{r, t) exp[i6(r, t)]. 

Since the definition l|12fl of the mean field involves a 
weighted spatial average over infinitely many local os- 
cillators, this quantity is expected to be smooth in space 
even these oscillators are behaving incoherently. This 
property of W is also clear from the differential form of 
Eq. Ull), i.e., 

0= -W + V 2 VF + exp(i0). (13) 

The above equation implies a strong similarity of W to 
B governed by Eq. J3J . If the mean field pattern rotates 
steadily with frequency Q, then R is time-independent 
and the relative mean field phase Oo defined by 

0(r,t) = nt + Q Q {r,t) 

is also time-independent. 

In terms of R(r,t) and <d(r,t), Eq. I|ll|) may be ex- 
pressed in the form of a one-oscillator dynamics 

d t cj){r, t)=uj- R(r, t) sin (<p(r, t) + a - 8(r, tj) , 

or if we introduce a relative phase variable tj)(r,t) 
through 

«£(r,t) = m + V(r,i), 

we have 

d t ip(r, t)=w-Cl- R{r, t) sin (ip(r, t) + a- 8 (r, t)) . 

(14) 



The definition of the mean field given by Eq. I|12|) be- 
comes 

R(r, t) exp[i© (r, t)] = J G(r - r') exp[iip(r' , t)]dr'. 

(15) 

Note that the set of Eqs. (|14fl and (|15|) is still equivalent 
to the original phase equation (|TT1) . 

We now proceed to some anatomy of the anomalous 
core structure taking advantage of the numerically ob- 
served fact that the mean field pattern has a well-defined 
center of rotation (chosen to be r = 0) at which W = 0. 
One may thus imagine a linear cross section S of the pat- 
tern passing through r = and study the radial profiles 
of various quantities emerging along S. Some results ob- 
tained in this way of analysis are summarized in FIG. [5] 
(a) to (c). 

An instantaneous radial profile of the mean field mod- 
ulus R is presented in FIG. El (a). As expected, it has a 
vanishing value at the origin, and its temporal fluctuation 
is also found negligibly small. 

Figure|5l (b) shows an instantaneous distribution of the 
phases <f> of the local oscillators lying on S (indicated by 
crosses ). The same panel also includes the pattern of 
the mean field phase on the same cross section (in- 
dicated by open circles). It is clear that there exists a 
well-defined critical radius separating the domains of co- 
herent and incoherent oscillators from each other. We 
also confirmed (but not shown explicitly) that the pro- 
files of the mean field phase and that of the phases of 
the coherent oscillators are almost stationary except for 
a vertical drift with constant velocity f2. 

Our interpretation of the results given by FIGs. (a) 
and (b) is that the entire system now splits into two 
subdomains such that the oscillators in one domain syn- 
chronize completely with the periodic mean-field forcing, 
while those in the other domain fail in synchronization. 
Further evidence supporting this interpretation is pro- 
vided in FIG. (c) where the distribution of the mean 
frequency Q (defined by a long-time average of dt<f>) of 
the local oscillators lying on S is shown. This frequency 
pattern is clearly composed of two parts. Namely, in the 
outer domain the oscillators have an identical frequency, 
while in the inner domain the frequencies are distributed, 
the latter implying phase randomization consistent with 
the scattered dots appearing in FIG.El(b). 

In the next section, we develop a theory for determin- 
ing the mean field pattern together with its rotation fre- 
quency, and also the motion of the individual oscillators 
driven by this mean field, in a self-consistent manner. 

IV. THEORY 

The basic equations to work with are Eqs. (|14fl and 
(|15fl . Our theory starts with the assumption that the 
mean-field pattern is steadily rotating, and therefore we 
drop the i-dependence from R and Oo in these equa- 
tions. A complete solution to this system of equations 
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FIG. 4: Spiral pattern (left) and its core structure (right) 
exhibited by non-locally coupled phase oscillators governed 
by Eq. ijTT). where a = 0.3. 

can be obtained in the following two steps. We first solve 
Eq. I|14|) for each ip as a function of R and Qq, which is 
easy to do. Note that R and 9o are the quantities yet 
to be determined. Secondly, the entire set of these so- 
lutions is substituted into Eq. (JTSJ) . The right-hand side 
of Eq. I|15|l thus becomes a functional of the mean field. 
In this way, the mean field value at each spatial point is 
expressed in terms of a functional of the mean field it- 
self. Solution of this functional self-consistency equation 
exists only for a special value of the rotation frequency f2 
of the mean field pattern. We will therefore be working 
with a non-linear eigenvalue problem. The final solution 
of this functional equation could be found only numeri- 
cally. 

The above self-consistent way of finding a solution 
to a many-oscillator problem resembles strongly Ku- 
ramoto's 1975 theory of synchronization phase transi- 
tion in a large population of globally coupled oscillators 
with distributed natural frequencies |9j- The main differ- 
ence is that the oscillators are now coupled non-locally 
rather than globally, and consequently the mean field is 
generally space-dependent leading to a functional self- 
consistency equation rather than a simple transcendental 
equation. Although the natural frequencies of the oscilla- 
tors are identical in the present case, the actual frequen- 
cies can be distributed due to the existence of a spatial 
gradient of the mean field. A simpler, one-dimensional 
version of the present type of theory based on a sim- 
ilar model of non-locally coupled phase oscillators was 
reported earlier pfj). 

An important feature common to all such theories is 
that the one-oscillator equation which involves the mean 
field amplitude as a parameter admits either a station- 
ary solution or a drifting solution. Which one to hold 
depends on the modulus of the mean-field. The crucial 
point to the theory is how to deal with the drifting so- 
lutions, because a simple substitution of this type of so- 
lutions into the definition of the mean field apparently 
contradicts the assumed stationarity of the mean field 
(in a suitable co-moving frame of reference). The seeming 
contradiction here can be resolved by using the invariant 
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FIG. 5: Radial profiles of various quantities corresponding to 
the spiral core of FIG. |1] (a) Instantaneous radial profile of 
the mean-field modulus R. (b) Instantaneous radial profile of 
the phases <fi of the local oscillators (crosses) and that of the 
mean field phase O (open circles), (c) Radial profile of the 
mean frequency O (defined by a long-time average of 9t0) of 
the local oscillators. 



measure associated with the drift motion. We will now 
show explicitly the steps leading to an exact solution to 
the problem. 

As stated above, there are two possible cases regarding 
the solution of Eq. 12}. They are: (Case I) \v — fi| < 
R, and (Case II) \w — fi| > R. Correspondingly, the 
oscillators are divided into two groups. In the first case, 
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which corresponds to the group of coherent oscillators, 
Eq. (|14f) admits a pair of stable and unstable fixed points. 
The stable one, denoted by ipo(r), is given by 



ip (r) = Sin 



R(r) 



e (r) 



The actual frequencies ui of the oscillators in this group 
are of course identical with fi. We substitute the above 
solution for ip(r) into Eq. 1)150. and restrict the integral 
to the domain where the inequality \uj — 0| < R(r) is 
satisfied. In this way, the contribution to the local mean 
field value coming from the coherent group of oscillators 
is obtained. 

The second case corresponds to the group of incoherent 
oscillators, for which Eq. (|14|) admits a drifting solution. 
The actual frequencies co{r) are now distributed and they 
are easily calculated as 



6j = n + 2n 



= n + (uj-n)\ i- 



R 



- o 



The contribution to the local mean field value from this 
incoherent group of oscillators can be found in the fol- 
lowing way. Since i\> is drifting, the factor exp(i?/0 in the 
integrand in Eq. (|15|l does not have a definite value. We 
are thus led to the idea that this factor should rather be 
replaced with its statistical average which can be calcu- 
lated by using the invariant measure, i.e., the probability 
density p(ip) associated with the drift motion. Noting 
that the probability density for the oscillator's phase to 
take on value ip must be inversely proportional to the 
drift velocity given by the right-hand side of Eq. (|14|) , we 
have 



p{ip) = C [uj - ft - i?sin (tp + a - 6 ) 



(16) 



where C is the normalization constant given by C — 

(2tt)^ 1 ( w - n)y/i - r?I{uj - n) 2 . 

Putting together the above-stated two types of contri- 
butions to the mean field, we finally obtain a functional 
self-consistency equation in the form 

R{r)e tQo{r) = J G(r- r')h(R(r'), 6 (r'), w - tydr', 

(17) 



where 



exp[iip (R,9 ,u- ^)] 



: ~ - ' ^j:^,R, e , w - n) ex P (* 

{\w-Cl\>R), 
or more explicitly 



' LU — il 

R 



x l-i 1- 



R 



lu — n 



R 



Numerical solution of Eq. (|17fl can be found iteratively. 
We did this in a finite domain defined by x,y € [0,40] 
with G appropriate for the free boundary conditions im- 
posed on Eq. (|13l) . Since a solution of Eq. (|17fl would 
only be available for a special value of Q — to which is still 
to be determined, its trial value was adjusted in each iter- 
ation step in such a way that a suitably defined distance 
between the two mean-field patterns, one produced at 
the current step and the other at the next step, may be 
minimized. In this way, by starting with a suitable ini- 
tial mean field pattern similar to the one obtained from 
numerical simulations, a rapid convergence of the mean 
field pattern and the value of was achieved. 

In FIG. our theoretical results obtained in this way 
are compared with the data given in FIG. 03 i.e., the 
results from direct numerical simulation of Eq. I|ll|l . The 
agreement is so excellent that our theory is expected to 
hold exactly in the continuum limit. 



V. SUMMARY AND CONCLUDING REMARKS 

Spontaneous generation of a local group of phase- 
randomized oscillators near the center of a rotating spi- 
ral pattern was confirmed through numerical simulations 
on non-locally coupled oscillators. It was argued that 
smaller value of the coupling strength favors the occur- 
rence of the core anomaly. The critical coupling strength 
K c associated with the onset of this anomaly was esti- 
mated from a simple argument. When K is sufficiently 
small, by which the oscillation amplitude even near the 
center of rotation takes almost a full value, a group of 
incoherent oscillators always exists. Still the overall spi- 
ral pattern looks completely normal. Guided by this 
fact observed numerically, we applied the phase reduction 
method for the purpose of gaining a deeper understand- 
ing of the phenomenon. The resulting phase oscillator 
model with non-local coupling was found to exhibit the 
same type of core anomaly. Under the assumption that 
the pattern of a suitably defined mean field is steadily ro- 
tating in spite of the existence of incoherence, we derived 
a functional self-consistency equation to be satisfied by 
the mean field. Its solution successfully reproduced var- 
ious results obtained from our direct numerical simula- 
tions carried out on this phase model. 

Finally, we remark that the present study is confined 
to a particular domain of parameter values where the 
mean field dynamics is regular. Our preliminary study 
suggests that under different conditions more complex 
collective dynamics occurs, which is characterized, e.g., 
by an elongation of the domain of incoherent oscillators 
and its irregular motion^. For the case of non-locally 
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coupled FitzHugh-Nagumo oscillators for small K, this 
occurs for larger b, i.e., when the symmetry of the local 
oscillator dynamics is lowered, although a perfect sym- 
metry (b = 0, or the van der Pol limit) is not neces- 
sary for the steady rotation of the mean field. For larger 
K (still below K c ), in contrast, regular dynamics of the 
mean field and circular shape of the domain of incoher- 
ent oscillators seem to persist against relatively strong 
asymmetry of the oscillator dynamics. These results will 
be reported elsewhere. 





FIG. 6: Comparison between the theory and numerical simu- 
lation. Theoretical results are indicated with solid lines in (a) 
and (c), and solid lines and scattered dots in (b). Numerical 
data, which are the same as those given in FIG. 03 are indi- 
cated with open circles and crosses, (a) Instantaneous radial 
profile of the mean-field modulus R. (b) Instantaneous radial 
profile of the mean- field phase O, and that of the phases <j> 
of the local oscillators, where the theoretically obtained scat- 
tered dots are the random numbers chosen from the proba- 
bility distribution given by Eq. I|16[l . (c) Radial profile of the 
mean frequency Co of the local oscillators. 
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